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Abstract 

We discuss a method that allows to compute the thermodynamic Casimir 
force at a given temperature in lattice models by performing a single Monte 
Carlo simulation. It is analogous to the one used by de Forcrand and Noth 
and de Forcrand, Lucini and Vettorazzo in the study of 't Hooft loops and 
the interface tension in SU(N) lattice gauge models in four dimensions. We 
test the method at the example of thin films in the XY universality class. In 
particular we simulate the improved two-component model on the simple 
cubic lattice. This allows us to compare with our previous study, where we 
have computed the Casimir force by numerically integrating energy densities 
over the inverse temperature. 
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1 Introduction 



In 1978 Fisher and de Gennes pj realized that when thermal fluctuations are re- 
stricted by a container a force acts on the walls of the container. Since this effect is 
similar to the Casimir effect, where the restriction of quantum fluctuations induces 
a force, it is called "thermodynamic" Casimir effect. Since thermal fluctuations 
only extend to large scales in the neighbourhood of a continuous phase transitions 
it is also called "critical" Casimir effect. Recently this effect has attracted much 
attention, since it could be verified for various experimental systems and quantita- 
tive predictions could be obtained from Monte Carlo simulations of spin models [2]. 
The neighbourhood of the critical point implies that the Casimir force is described 
by a universal finite size scaling function 0. For the film geometry, finite size scaling 
predicts 

Fca~(i:o,t) ^ ^e{t[L^/ioY'^) (1) 

where /c^ is the Boltzmann constant, T the temperature, t = {T —T^) /T^. the reduced 
temperature and Lq the thickness of the film. The amplitude of the correlation 
length in the high temperature phase is defined by 

e ^ eot"^ (2) 

where ^ is the correlation length of the bulk system and v its critical exponent. The 
function 9{x) is the same for all films in a given universality class, where also the 
boundary universality class [9] has to be taken into account. 

As a first application of the numerical method discussed here, we study the 
improved two-component 0^ model on the simple cubic lattice. The phase transition 
of this model belongs to the XY universality class in three dimensions. Also the 
A-transition of ^He shares this universality class. The experimental study of the A- 
transition provided highly accurate estimates for critical exponents and amplitude 
ratios of the bulk system. For a review see [10]. Also confined systems have been 
studied in detail at the A-transition of "^He [11]. In particular, the thermodynamic 
Casimir force in thin films of "^He has been measured [T^ [T3] . These experiments 
confirm that the thermodynamic Casimir force for films of different thickness Lq 
can indeed be described by the same scaling function 9{x). For all temperatures the 
force turns out to be negative. In the high temperature phase 9{x) is monotonically 
decreasing with decreasing x. The Casimir force vanishes for large values of x. 

"'^For reviews on finite size scaling see [31 For reviews on critical phenomena and the Renor- 
malization group see e.g. [3 [6l [71 [8] . 
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At the critical point of the bulk system 9(0) = -0.07 ± 0.03 [T2]. In the low 
temperature phase the finite size scaling function shows a minimum at Xmin ~ —5.5 
with Ojjiin ~ 1-3 [H]. For x < Xmin the finite size scaling function increases with 
decreasing temperature. For small values of x it seems to approach a finite negative 
value. 

It has been a challenge for theorists to compute the finite size scaling func- 
tion 9{x). Krech and Dietrich p3l IT5] have computed it in the high temperature 
phase using the e-expansion up to 0(e). This result is indeed consistent with the 
measurements on ^He films. Deep in the low temperature phase, the spin wave 
approximation should provide an exact result. It predicts a negative non- vanishing 
value for 9{x). However the experiments suggest a much larger absolute value for 
9{x) in this region. Until recently a reliable theoretical prediction for the minimum 
of 9{x) and its neighbourhood was missing. Using a renormalized mean-field ap- 
proach the authors of [121 [IZ] have computed 9{x) for the whole temperature range. 
Qualitatively they reproduce the features of the experimental result. However the 
position of the minimum is by almost a factor of 2 different from the experimental 
one. The value at the minimum is wrongly estimated by a factor of about 5. 

Only quite recently Monte Carlo simulations of the XY model on the simple 
cubic lattice [18], UHl [20] provided results for 9{x) which essentially reproduce the 
experiments on ^He films [121 US] • In [21] we have applied the method used by [I9] 
to study the improved two-component 0^ model on the simple cubic lattice. The 
study of this model should provide more accurate results since corrections oc Lq 
with uj = 0.785(20) [22] are eliminated. Essentially our result confirms those of 
[T8| [HI [2^. However there is a discrepancy in the position Xmin of the minimum of 
9{x) that is clearly larger than the errors that are quoted: In Xmin = —5.3(1) and 
in [20] Xmin = —5.43(2) which has to be compared with our result Xmin = —4.95(3) 

In order to verify our result [21] we compute the thermodynamic Casimir force 
in the two-component (p^ model using a different method that is analog to that of 
[221 [21] used to compute the string tension and 't Hooft loops in lattice gauge model. 
The general idea is similar to that of [TSl [20]- However, in contrast to [TSl [20] a 
single simulation is sufficient to obtain the Casimir force at a given temperature. 

This paper is organized as follows: First we define the 0^ model on the simple 
cubic lattice. Then we discuss in detail the method used here to compute the ther- 
modynamic Casimir force. In section [H we discuss our numerical simulations. First 
we performed numerical simulations at the critical point of the three-dimensional 

^Provided that fbuik is known. 
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system. Next we computed the free energy density for the thermodynamic hmit of 
the three-dimensional system at two values of the inverse temperature j3 in the high 
and the low temperature phase each. Then we have measured the thermodynamic 
Casimir force for Lq = 8.5 at various temperatures. Finally we have simulated at 
Xmin for the thicknesses Lq = 6.5, 7.5, 9.5, 12.5 and 24.5 to complement our results 
of ref. [21]. Finally we summarize our results and give our conclusion. 



2 The model and the observables 

We study the two-component (j)'^ model on the simple cubic lattice. We label the 
sites of the lattice by x = (xq, xi, X2). The components of x might assume the values 
Xi G {1,2,..., Li}. We simulate lattices of the size Li = L2 = L and Lq <^ L. In 
1 and 2-direction periodic boundary conditions are used. In order to mimic the 
vanishing order parameter that is observed at the boundaries of ^He films, free 
boundary conditions in 0-direction are employed. This means that the sites with 
Xq = 1 and Xq = Lq have only five nearest neighbours. This type of boundary 
conditions could be interpreted as Dirichlet boundary conditions with as value of 
the field at = and xq = Lq + 1. Note that viewed this way, the thickness of 
the film is Lq + 1 rather than Lq. This provides a natural explanation of the result 
Ls = 1.02(7) obtained in [25]. The Hamiltonian of the two-component 0^ model, 
for a vanishing external field, is given by 

n = -/3j2 ^^■$y + T.[^' + A(0^ - 1)2] (3) 

<x,y> X 

where the field variable (p^ is a vector with two real components. < x,y > denotes 
a pair of nearest neighbour sites on the lattice. The partition function is given by 



n 



1) / ^^(2) 



exp{-n). (4) 



Note that following the conventions of our previous work, e.g. [26] , we have absorbed 
the inverse temperature /3 into the Hamiltonian. [f] In the limit X ^ 00 the field 
variables are fixed to unit length; Hence the XY model is recovered. For A = we get 
the exactly solvable Gaussian model. For < A < cxd the model undergoes a second 
order phase transition that belongs to the XY universality class. Numerically, using 



^Therefore, following J] we actually should call it reduced Hamiltonian. 
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Monte Carlo simulations and high-temperature series expansions, it has been shown 
that there is a value A* > 0, where leading corrections to scaling vanish. Numerical 
estimates of A* given in the literature are A* = 2.10(6) [27], A* = 2.07(5) [26] and 
most recently A* = 2.15(5) |22j|. The inverse of the critical temperature I3c has been 
determined accurately for several values of A using finite size scaling (FSS) [22] . 

We shall perform our simulations at A = 2.1, since for this value of A compre- 
hensive Monte Carlo studies of the three-dimensional system in the low and the 
high temperature phase have been performed [251 [221 1211 129]. At A = 2.1 one gets 
(3c = 0.5091503(6) [22]. Since A = 2.1 is not exactly equal to A*, there are still cor- 
rections oc Lq'^, although with a small amplitude. In fact, following [22], it should 
be by at least a factor 20 smaller than for the standard XY model. 

2.1 The energy density and the reduced free energy 

Note that in eq. ([3]) (3 does not multiply the second term. Therefore, strictly 
speaking, (3 is not the inverse of ksT. In order to study universal quantities it is 
not crucial how the transition line in the jS-X plane is crossed, as long as this path is 
not tangent to the transition line. Therefore, following computational convenience, 
we vary /3 at fixed A. In the following equations it is understood that A is kept 
fixed. 

The reduced free energy density is defined as 



Note that compared with the free energy density /, a factor ksT is skipped. For 
convenience we have defined the reduced free energy such that /(O) = 0. For j3 = 
the partition function factorizes and thus log Z {0) / {LQL1L2) does not depend on 
the system size. 

We define the (internal) energy density as the derivative of the reduced free 
energy density with respect to f3. Furthermore, to be consistent with our previous 
work e.g. [30], we multiply by —1: 



1 



[logZ(/3)-logZ(0)] . 



(5) 



E 



1 dlogZ 
L0L1L2 d(3 



(6) 



It follows 
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which can be easily determined in Monte Carlo simulations. From eqs. (I5ll6!) it 
follows that the free energy density can be computed as 

f{P) = fiPo)- ^dpEiP) . (8) 

3 The numerical method 

From a thermodynamic point of view, the Casimir force per unit area is given by 

Fcasimir kgT (9) 

where Lq is the thickness of the film and f^x = ffUm — Lofbuik is the reduced excess 
free energy per area of the film. In lattice models the thickness Lq assumes only 
integer values. Therefore we have to approximate the derivative by a finite difference 
Fcasimir (Lo,t) ^ -kBTAfex{Lo,t), where 

A/e.(Lo, t) = /(Lo + 1/2, t) - f{Lo - 1/2, t) - h^ikit) (10) 

where Lq + 1/2 is integer. /(Lq + 1/2, t) and /(Lq — 1/2,/!:) are the reduced free 
energies of films of the thicknesses Lq + 1/2 and Lq — 1/2, respectively, and fbuikii) 
is the reduced free energy density of the three-dimensional bulk system. The main 
numerical task is to compute the difference of the reduced free energy per area for 
films of the thickness Lq — 1/2 and Lq + 1/2. 

In order to compute this difference, it is useful to have the same number of field- 
variables for both systems. To this end, we add Li x L2 isolated sites to the film of 
the thickness Lq — 1/2. Isolated means that the field at such a site is subject to 
the potential 0^ + A(0^ "1)^ but the interaction with other sites is missing. Using 
our definition ([5]), adding isolated sites to the film does not change the free energy 
per area. Let us denote the partition function of this system by ZiQ_i/2. Now we 
can express the difference of the reduced free energies as 

F(Lo + l/2,t) - F(Lo-l/2,t) = log|^^ 

D[0] exp(-ffL,+i/2) exp(-/? E<x,y>G[Lo+i/2] '^^■'^y) 
= log 



D[0] exp(-/7L„+i/2) 

log(exp(-/? $x$y)) , (11) 

Lo+1/2 
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where < x,y >G [Lq + 1/2] denotes the sum over all nearest neighbour pairs, 
where at least one of the sites is an element of the layer xq = + 1/2. Formally 
we have rewritten the difference of free energies as an expectation value. The 
problem is that the observable is strongly fluctuating and therefore it is impossible 
to obtain an accurate estimate from a Monte Carlo simulation of the film of the 
thickness Lq + 1/2. A well known method to overcome this problem is the so 
called "multistage sampling" strategy; see e.g. [31j. This means that a sequence of 
systems is introduced that interpolates between the two we are interested in. These 
systems are characterized by the Hamiltonians Hq, Hi, . . . , if at, where we identify 
Hq = Hlo-i/2 and Hn = -f^Lo+i/2- Now we can rewrite the ratio of partition 
functions as 

Zn Zi Z2 Zn 
where we can write the factors as 

^. 

Zi+i = — ^ = (exp(-ifi+i + Hi)),+i (13) 

and hence 

N 

F{Lo + 1/2, t) - F{Lo - 1/2, t) = J2 log ■ (14) 

i=l 

If the sequence is properly chosen and N sufficiently large, the fluctuations of 
exp(— + Hi) are small and the expectation value can be accurately determined 
from the simulation of the system z + 1. @ Obviously there is much freedom in the 
construction of the sequence of systems. A straight forward one is given by 

H, = Hl,-1/2 + ^P J] $Jy. (15) 

<x,s/>e[Lo+l/2] 

This choice is very similar to the one used by [20] . The main difference is that 
these authors did consider as starting system a film of thickness Lq — 1/2 plus a 
two-dimensional system of the size Li x L2. This means that in contrast to our 
choice the intra-layer couplings are switched on. 

Here we use a different interpolation. It is inspired by a method used to compute 
the string tension and 't Hooft loops in lattice gauge theories [23], [21] . 

■*It might be even better to express the difference as an expectation value in a system that is 
between i and i + 
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Figure 1: We sketch the layer = I/q + 1/2 of our system. In the sketch Li = 
L2 = 6. The sites are given by circles. The filled ones are coupled to the system, 
while the empty ones are isolated. We compute the free energy difference between 
the system, where the shaded circle is isolated and the system, where it is coupled 
to the film. For a discussion see the text. 



We add the isolated sites one by one to the film. In the step i = X1L2 + X2 the 
site X = {Lq + 1/2, xi, X2) is added. All sites that are added are coupled with their 
nearest neighbours that are already in the film. This way we have constructed a 
sequence of Li XL2+I systems. Hence, L1XL2 independent Monte Carlo simulations 
have to be performed to obtain F(Lo + 1/2, t) - F(Lo - 1/2, t). Following [231 El] 
this can however be avoided: With increasing Li,L2 the sum ( 1141) is dominated 
by contributions where the newly added site is far from the boundaries in 1 and 
2-direction. Hence most of the contributions are essentially equal to that for xi = 
Li/2 and X2 = L2/2 as sketched in figure [H In the limit Li,L2 — > 00, this should 
become exact. Hence only a single simulation for xi = Li/2 and X2 = L2/2 is 
required. 

Usually updates are performed on the whole lattice El before a measurement 
^In the Monte Carlo slang this is often called "sweep" . 
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of the observables is performed. However in the present case, the observable is 
locahzed at a single site. Therefore the effort for the measurement and the update 
would be highly unbalanced. To circumvent this problem, one would like to update 
the fields in the neighbourhood of this site more frequently than those far off. In 
order to achieve this we follow the idea presented in [32]: We consider a sequence 
of sub-sets of the sites of the lattice. Here, the smallest set consists of the site 
(Lq + 1/2, Li/2, L2/2) only. The next larger one consists of (Lq + 1/2, Li/2, L2/2) 
and its three neighbours (Lq + l/2,Li/2 - 1,^2/2), (Lq + 1/2, Li/2, L2/2 - 1) 
and (Lq — 1/2, Li/2, L2/2). The larger ones are given by blocks of the size bi x 
(26/ + 1) X {2bi + 1) and Lq x (26/ + 1) x (26/ + 1) if 6/ > Lq. These blocks are 
centred around the site {Lq + 1/2, Li/2, L2/2). If 6/ < Lq, the eight corners of 
these blocks are {Lq + 1/2, Li/2 - 6/, L2/2 - 6/), {Lq + 1/2, Li/2 - 6/, L2/2 + 6/), 
(Lo-l/2,Li/2-6/,L2/2 + 6/), (Lq - 1/2, Li/2 + 6/, L2/2 + 6/), (Lq + 1/2 - 6/, Li/2 - 
6/, L2/2 - 6/), {Lq + 1/2 - 6/, Li/2 - 6/, L2/2 + 6/) , {Lq + 1/2 - 6/, L,/2 - 6/, L2/2 + 6/) 
and {Lq + 1/2 — 6/,Li/2 + 6/,L2/2 + 6/). In our simulations we have used 6/ = 
1, 2, 3, 5, 10, 20, 40, 80, . . ., where the largest 6/ is chosen such that 6/ < Li, L2. 

In a certain sequence. Metropolis and overrelaxation sweeps |§ over these sub- 
sets are performed. This sequence, which we shall call one update cycle, is best 
explained by the following pseudo-code 

cluster_update ; metrosweep(full lattice); oversweep(f ull lattice); 
for(il=0; il < m.l; il++) 
{ 

metrosweep(b_l) ; oversweep(b_l) ; 
for(i2=0; i2 < m_2; i2++) 
{ 

metrosweep(b_2) ; oversweep(b_2) ; 



for(iM=0; iM < m_M; iM++) 
{ 

metrosweep(b_M) ; oversweep(b_M) ; 

measure ; 

} 



^We have implemented these updates as discussed in appendix A of [26] 
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} 

} 



We did not perform an accurate tuning of the parameters mi,m2,m^, . . .niM] 
instead we have chosen them such that the CPU-time spent at each block-size is 
roughly the same. 

In the case of the single cluster-updates ^33j it is easy to focus on the site (Lq -|- 
1/2, Li/2, L2/2). One simply starts the clusters at the site (Lq -|- 1/2, Li/2, L2/2) 
instead of choosing the starting point at random. In our numerical tests we have 
not yet implemented this idea. 

3.1 The measurement 

The measurement consists in its most naive implementation in the evaluation of 

A = exp(-/30(io+i/2,Li/2,L2/2) ■ *&(Lo+l/2,Li/2,L2/2)) (16) 

where 

'^'(Lo+l/2,Li/2,L2/2) = 0(Lo+l/2,Li/2-l,L2/2) + 0(Lo+l/2,Li/2,L2/2-l) + 0(Lo-l/2,Li/2,L2/2) 

(17) 

We have reduced the variance by performing the integral over the angle of the field 
0{Lo+i/2,Li/2,L2/2) exactly. This results in the improved observable 



/g'^da exp(-i2cosa) 27rJo(-R) 

where 

R = /5|0(Lo+l/2,Li/2,L2/2)||*(Lo+l/2,Li/2,L2/2)| (19) 

and /o(-R) is a modified Bessel function. For our simulations we have tabulated 
l/(27r/o(i?)) for < i? < 3 with a step-size of 0.0001, i.e. for 30001 values of R. 
During the simulation we then evaluated l/(27r/o(-R)) for < i? < 3 by quadrat- 
ically interpolating the results given in the table. If i? > 3 we have evaluated the 
integral in eq. f|T8l) numerically. A sufficient precision can already be achieved with 
about 30 nodes. 

The expectation value 

z={A) = {A) (20) 
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is estimated by averaging A over all measurements that we performed after tlier- 
malization. Here we have dropped the subscript i = (Li/2)L2 + L2/2 of eq. [131 
since only this value of i will be considered in the following. During the simulation 
we have averaged already all measurements in a given update cycle. These averages 
were written to a file. The statistical error was then computed taking into account 
the integrated autocorrelation time of these cycle averages. 

4 Numerical Results 

First we have simulated at the critical temperature of the bulk system. Next we have 
determined the reduced free energy of the bulk system at /5 = 0.49 and (3 = 0.5 in 
the high temperature phase and at /3 = 0.533 and (3 = 0.56 in the low temperature 
phase. For Lo,Li,L2 ^ C, the reduced free energy of the bulk system is given by 
fbuik = log z. Our new results are consistent with those obtained by integrating 
the energy densities computed in [50]. Then we have studied films of the thickness 
Lq = 8.5 at four temperatures in the low temperature phase of the bulk system. 
Also here we find that the results are consistent with those of [30]. Finally we 
have simulated the thicknesses Lq = 6.5, 7.5, 9.5, 12.5 and 24.5 at Xmin- These 
simulations complement our results of [30] at Xmin- 

As random number generator we have used the SIMD-oriented Fast Mersenne 
Twister algorithm [31]. 

4.1 Simulations at the critical point 

First we performed simulations at the inverse critical temperature (3c = 0.5091503(6) 
of the three-dimensional system using lattices of the thicknesses Lq = 8.5, 12.5, 
16.5, 24.5, 32.5 and 64.5. In all cases we have chosen Li = L2 = 12.5 x {Lq — 
1/2). Since the correlation length of the film is ^2ndjUm/Lo,eff ~ 0.416 [30] this 
should be sufficient to keep deviations from the two-dimensional thermodynamic 
limit smaller than our statistical errors. As a check we have simulated for Lq = 8.5 
in addition Li = L2 = 20,30 and 50. We find z = 0.84950517(36), 0.84951362(36) 
and 0.84951572(37) for these lattice sizes, respectively. Indeed, starting from Li = 
L2 = 30 our results are consistent within error bars. Our results for Li = L2 = 
12.5 X [Lq — 1/2) are summarized in table [H In these simulations we have used 
block sizes up to bi = 10, 20, 20, 40, 40 and 80 for Lq = 8.5, 12.5, 16.5, 24.5, 32.5 
and 64.5, respectively. For all these thicknesses and for all block-sizes we have used 
nil = 6. The number of update cycles is 21349600, 7738500, 7700000, 2040000, 
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Table 1: Results for z at [3c = 0.5091503 (22] for lattices of the size Li = L2 = 
12.5 X Lq. 





z 


8.5 


0.84951552(24) 


12.5 


0. 


84947657(16) 


16.5 


0. 


84946525(23) 


24.5 


0. 


84945897(13) 


32.5 


0. 


84945717(16) 


64.5 


0. 


84945602(14) 



1308200 and 281700 for Lq = 8.5, 12.5, 16.5, 24.5, 32.5 and 64.5, respectively. In 
total these simulations took about 16 month of CPU-time on a single core of a 
Quad-Core Opteron(tm) 2378 CPU (2.4 GHz). 
The reduced excess free energy behaves as 

fULo, t) = L-ljfh{t[L,li,Y/^) + hit) (21) 

where Lo^^ff = Lq + Lg with Lg = 1.02(7) ^23] takes into account corrections due to 
the Dirichlet boundary conditions and b{t) gives the effect of the Dirichlet boundary 
conditions on the analytic part of the free energy of the film. For a discussion and 
references see [35]. Taking the derivative with respect to Lq at t = we arrive at 

= 2^(0)^01// = ^(0)^0";'// (22) 

t=o 

where 6 is the finite size scaling function of the thermodynamic Casimir force. 
It follows 

iog^(Lo,/5j = - om^.iff. (23) 

Note that in the thermodynamic limit the singular part of the free energy density 
vanishes at the critical point; hence fbuikiPc) = fnsiPc)- The results of our fits are 
given in table O 

In order to estimate the effect of the error of Ls on our results we have repeated 
these fits using Lg = 0.95. E.g. for Lo^rnin = 12.5 we get fns = —0.16315930(10) 
and ^(0) = —0.0593(5). We have also checked the effect of the error of f3c- To this 
end we have computed A/(Lo, 0.5091509) by using the data for the energy given in 



dfex{Lo,t) 
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Table 2: Results for Gts with the ansatz / f23)) . where we have used Lg = 1.02 as 
input. All data for Lq > Lo^min are fitted For a discussion see the text. 



f-'Q,inin 


fns 


^^(0) 


xVd.o.f. 


8.5 


-0.16315935(9) 


-0.0606(3) 


0.20 


12.5 


-0.16315932(10) 


-0.0603(6) 


0.14 


16.5 


-0.16315930(12) 


-0.0597(17) 


0.15 



table 1 of We find that the effect on f^s and 9{0) is small and can be ignored 
here. Based on the result obtained for Lo^min = 12.5 we take as final results 

fns = -0.1631593(1) , 6(0) = -0.060(2) (24) 

where the error-bar covers both the statistical error as well as the error due to the 
uncertainty of L^. 

This can be compared with the result for ^He films ^(0) = —0.07 ±0.03 [12], the 
e-expansion up to 0(e): ^^(0) = —0.044 taken from table I of [2] and the estimate 
^(0) = —0.062(5) obtained from Monte Carlo simulations of the standard XY model 
[T9] . The authors of [36] quote h{0) ^ —0.03 (in their notation A-^) as final result. 
All these results are consistent with ours. The largest discrepancy is seen for the 
e-expansion. However one should note that in |21] it has been observed that 
in the high temperature phase for x ^ 1 the numerical result for 6 matches nicely 
with the e-expansion [14j . 

4.2 Free energy density of the bulk system 

Here we compute the free energy density of the bulk system for two values of (3 
in the high temperature phase and two values of (3 in the low temperature phase. 
These results are compared with ones obtained by integrating the energy density 
starting from Pc = 0.5091503 using the start value fiPc) = —0.1631593(1) obtained 
above. 

For sufficiently large Lq, Li and L2 the quantity log z should be a good approxi- 
mation of the bulk free energy density. In particular in the high temperature phase, 
this should be the case for Lo,Li,L2 3> ^3d- Here we performed simulations at 
p = 0.49 where 6nd,3Z) = 3.72370(19) and p = 0.5 where 6nd,3D = 6.1498(5) (see 
table 5 of 
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At /3 = 0.49 we have simulated Lq = 49.5, Li = L2 = 50 and Lq = 99.5, 
Li = L2 = 100. For Lq = 49.5, Li = L2 = 50 we have used block sizes up to 
bi = 20 and mi = 6. From 5203000 cycles we get /(0.49) = -0.14712079(18). For 
Lq = 99.5, Li = L2 = 100 we have used block sizes up to bi = 40 and mi = 6. 
From 945000 cycles we get /(0.49) = -0.14712095(17). As expected, these resuhs 
are indeed consistent within error bars and hence a good approximation of the 
thermodynamic limit. 

Based on the experience gained at = 0.49 we have simulated ai (3 = 0.5 only 
the lattice size Lq = 99.5, Li = L2 = 100. We have used block sizes up to 61 = 20 
and mi = 6. From 598000 cycles we get /(0.5) = -0.15519942(24). 

In the low temperature phase we find from simulations of a 199.5 x 500^ lattice 
/(0.533) = -0.18931867(66) and /(0.56) = -0.22693625(73). We have used blocks 
up to the size bi = 80 and mi = 6 for all block sizes. We performed 24700 and 
21700 cycles for j3 = 0.533 and P = 0.56, respectively. Both of these simulations 
took about 8 weeks of CPU-time on a single core of a Quad-Core Opteron(tm) 2378 
CPU (2.4 GHz). 

Now we can check whether these results for the free energy density are consistent 
with those obtained from integrating the energy density |30] using eq. ([3]). 

In [30] we have computed the energy density of the three-dimensional bulk sys- 
tem in the range of inverse temperatures 0.49 < (3 < 0.58. We have fitted these 
data in the range 0.49 < (3 < 0.529 with the ansatz 

E{(3) = Ens + Cns{(3 - (3c) + a±\(3- + dns{(3 - (3cf + b±\(3- (25) 

where Ens, Cns, Pc = 0.5091503(6) and a = —0.0151(3) [22] are input and a±, dns 
and b± are the 5 free parameters of the fit. For (3 < 0.529 we have integrated this 
ansatz, using the results for the fit-parameter obtained in [30] . In all cases we have 
taken Po = Pc = 0.5091503 as starting point of the integration, where we have used 
the estimate of f{Pc) obtained above. Our results are summarized in table [31 For 
P > 0.529 we performed a numerical integration of the energy density using the 
trapezoidal rule, starting from Po = 0.52. The estimate for /(0.52) is taken from 
table O We have checked that our result virtually does not depend on the choice 
of Po, where we switch from the integration of the ansatz fl25|) to the numerical 
integration of the energy density. Also the results for P > 0.529 are given in table 
ini The error quoted is dominated by the error for the free energy at Pc- 

In table [3] we also give our results for the free energy density of the bulk system 
at P = 0.49, 0.50, 0.533 and 0.56 as computed by the new method discussed here. 
We find that the results are consistent within error-bars. This confirms that we 
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Table 3: Numerical results for the free energy density of the bulk system. These 
where obtained by integration of the energy density. As starting point of the in- 
tegration we have taken the critical point j3c and the value f{f3c) obtained in the 
previous subsection. In addition in the third column we give estimates of the free 
energy density obtained directly with the method discussed in the present work. 



p 


/ INTEGRAL 


/ DIRECT 


0.49 


-0.1471210(1) 


-0.1471210(2) 


0.50 


-0.1551994(1) 


-0.1551994(2) 


0.51463 


-0.1684460(1) 




0.52 


-0.1740847(1) 




0.52348 


-0.1779533(1) 




0.5301 


-0.1857405(1) 




0.533 


-0.1893183(1) 


-0.1893187(7) 


0.53814 


-0.1958961(2) 




0.54 


-0.1983478(2) 




0.54432 


-0.2041851(2) 




0.56 


-0.2269362(2) 


-0.2269363(7) 



can indeed compute the free energy density of the bulk system with 6 to 7 accurate 
digits. 

4.3 Films of the thickness Lq = 8.5 

We have simulated at = 0.52, 0.533, 0.54, 0.56 in the low temperature phase of the 
three-dimensional system. We have taken lattices of the size Li = L2 = 50, 100, 250, 
500 and 1000 to control corrections to the two-dimensional thermodynamic limit 
of the thin film. The simulations for L = 1000 took about 18 days of CPU time 
each. In table H] we give our results for —Afex, where we have used the estimates 
of fbuik computed in section 14.21 At P = 0.52 the results for all choices of Li = L2 
are consistent within error-bars. At f3 = 0.533 a clear deviation of the results from 
those obtained for larger lattices can be observed up to Li = L2 = 100. The result 
for L = 500 deviates by a bit more than two standard deviations from that for 
L = 1000, while the results for L = 250 and L = 1000 are consistent within error 
bars. For (3 = 0.54 and /3 = 0.56 the results obtained for L = 250, 500 and 1000 
are consistent within error-bars. We conclude that in all cases for L = 1000 the 
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Table 4: Numerical results for minus the derivative of the excess free energy — A/ex 
of films of the thickness Lq = 8.5. In the last row we give the results of I21]j for 
comparison. 



L\(3 


0.52 


0.533 


0.54 


0.56 


50 
100 
250 
500 
1000 


-0.0007423(14) 

-0.0007432(8) 

-0.0007436(5) 

-0.0007439(3) 

-0.0007433(3) 


-0.0014878(23) 

-0.0015797(13) 

-0.0015845(8) 

-0.0015863(5) 

-0.0015846(5) 


-0.0010417(27) 
-0.0011867(18) 
-0.0012564(14) 
-0.0012679(10) 
-0.0012666(11) 


-0.0003621(23) 

-0.0003870(13) 

-0.0003934(8) 

-0.0003940(5) 

-0.0003945(6) 


ref. [21] 


-0.0007392(18) 


-0.0015795(24) 


-0.0012600(26) 


-0.0003874(28) 



deviation from the thermodynamic hmit is smaller than the error bar. 

For comparison we give in the last row results [21] which were obtained by 
numerical integration of Monte Carlo data for AE^x- We see that the results of [21] 
are about 0.000004 larger than our present ones. This deviation is about twice the 
statistical error. In ^1] we have started the integration at /5 = 0.49 for Lq = 8.5, 
setting A/ex(0.49) = 0. From the e-expansion [H] we get 9{x) ~ —0.0039 for 
X = t[Lo,e///^o]^^'^ corresponding to (3 = 0.49 and Lq = 8.5. Hence —A/ex = 
OLqIjtj ^ —0.0000041 which fully explains the difference observed in table HI 

4.4 The minimum of the Casimir force 

In [21] we have determined the position of the minimum of 6 for a large number of 
thicknesses of the film. To this end we have determined the zero of 

AEex(Lo, P) = E{Lo + 1/2, P) - E{Lo - 1/2, P) - EbuikiP) (26) 

where E{Lq + 1/2, /5) is the energy per area of a film of the thickness Lq + 1/2 
and Ebuik^P) the energy density of the three-dimensional bulk system. We had 
simulated at a few values of /? in the neighbourhood of Pmin- To get a prehminary 
estimate of Pmin we used the information gained already from the simulations for 
Lq = 8.5, 16.5 and 32.5, where we have simulated a large range of /3- values and the 
ansatz Pmin{Lo) — Pc — -^o,e//- These results are given in table O which we have 
copied from table 2 of [21]. In the present work, we have added the values of 9min 
for Lq = 6.5, 7.5, 9.5, 12.5 and 24.5 that were missing in [21]. To this end, we have 
simulated lattices of the size Li = L2 = 500 for Lq = 6.5 and 7.5, Li = L2 = 1000 
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Table 5: The position (3min of the minimum of the Casimir force and its value 
—^fex,min 9^ sl function of the thickness Lq. In the present work we have completed 
the table by adding Afex,min for Lq = 6.5, 7.5, 9.5, 12.5, and 24.5. 







A f ex, min 


6.5 


0.54432(2) 


-0.0032744(13) 


7.5 


0.53814(2) 


-0.0022305(11) 


8.5 


0.53354(2) 


-0.001582(3) 


9.5 


0.53010(2) 


-0.0011714(8) 


12.5 


0.52348(2) 


-0.0005468(6) 


16.5 


0.51886(2) 


-0.0002494(11) 


24.5 


0.51463(2) 


-0.0000803(3) 


32.5 


0.51279(2) 


-0.0000348(5) 



for Lo = 9.5 and Lq = 12.5, and Li = L2 = 2000 for Lq = 24.5. From these 
simulation we get log 2, while fhuik is taken from table [31 Our results for —Af^x are 
given in table O 

Let us briefly discuss the simulation of the L = 24.5 film: The simulations took 
about 2 month of CPU-time on a single core of a Quad-Core Opteron(tm) 2378 
CPU (2.4 GHz). We performed 33000 update cycles. We have used block sizes up 
to bi = 160 and mi = 6 for all block sizes. 

First we have fitted the results for — A/ex.mm given in the third column of table 
Elwith the ansatz 

Afex,min ^mm(-^0 ~l~ -^s) (2''') 

where Omin and Lg are the free parameters of the fit. Our results are summarized 
in table 

The x^/d.o.f. is smaller than 1 starting from Lo,mm = 8.5, where all data with 
Lq > Lo^min are included into the fit. We find Lg ^ 0.89 which is a bit smaller than 
our previous result Lg = 1.02(7) [25]. Note that already in [21] we observed that 
Lg = 0.95 apparently leads to a better matching of the data than Lg = 1.02. 

To check the possible effect of sub-leading corrections we have fitted our data 
also with the ansatz 

~ Afex,min = ^mini^ cLq'^)(Lq + Lg) ^. (28) 

Note that there a number of different corrections with an correction exponent close 
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Table 6: We have Rtted the minimum of the thermodynamic Casimir force with the 
ansatz ^27\] 







Ls 


xVd.o.f. 


6.5 


-1.299(2) 


O.J 


>49(5) 


2.64 


7.5 


-1.305(3) 


O.J 


>64(7) 


1.64 


8.5 


-1.313(5) 


O.J 


>89(13) 


0.89 


9.5 


-1.310(5) 


O.J 


>80(15) 


0.34 


12.5 


-1.312(9) 


O.J 


,88(33) 


0.50 



Table 7: We have fitted the minimum of the Casimir force with the ansatz ( f25)) 



Lo,min 




Ls 


c 


xVd.o.f. 


6.5 


-1.322(8) 


0.953(3) 


1.08(35) 


1.13 


7.5 


-1.320(10) 


0.945(5) 


0.97(61) 


1.40 



to 2. E.g. oc Lq^ with u' = 1.8(2) [37J or the restoration of the symmetries that 
are broken by the lattice. Our results are summarized in table [281 

Now the value Lg ~ 0.95 is fully consistent with our previous result [25j. As final 
result we quote 9min = —1-31(2), where we have estimated the systematic error by 
the difference of the two fits fl27|28p . This result fully confirms our previous estimate 
= -1.31(3) p. 

5 Summary and Conclusion 

We have discussed a method to compute the thermodynamic Casimir force in lattice 
models which is closely related with the one used by de Forcrand and Noth [22] and 
de Forcrand, Lucini and Vettorazzo [23] in the study of 't Hooft loops and the 
interface tension in SU(N) lattice gauge models in four dimensions. 

We have tested the method at the example of thin films of the improved two- 
component 0^ model on the simple cubic lattice. This model shares the XY uni- 
versality class with the A-transition of '^He. Therefore the Casimir force that is 
measured for thin films of ^He [T2| [13] should be governed by the same univer- 
sal finite size scaling function 6 as that computed from lattice models in the XY 
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universality class. 

Only quite recently 9 has been obtained from Monte Carlo simulations of the 
standard XY model on the simple cubic lattice [18], [HI [20]. This result is of par- 
ticular interest, since other theoretical methods do not provide us with accurate 
results for 9 for the whole range of the scaling variable x = t[Lo/^o]^'''^- Overall one 
finds a reasonable match between the experimental and Monte Carlo results. In [21] 
we have redone the Monte Carlo simulations using the improved two-component 0^ 
model on the lattice. It turns out that there is a discrepancy in the position Xmin 
of the minimum of 9{x): Xmin = —5.3(1) [19] and Xmin = —5.43(2) [JD] have to be 
compared with our result Xmin = —4.95(3) [2T] . 

The purpose of the present work is twofold: First we like to figure out the 
performance of the method and secondly we like to check and to complement the 
results of [21]. In particular: 

We have accurately computed the finite size scaling function of the thermody- 
namic Casimir force 9{0) at the critical point of the three-dimensional bulk system. 
Our result is consistent with the experimental result for ^He films [12] and previous 
Monte Carlo simulations [191 ISB] . On the other hand there is a clear discrepancy 
with the e-expansion [14]. 

We have demonstrated that the method even allows to compute the free energy 
density of the bulk system. However it seems to be more efficient in this case to 
integrate the energy density ([3]). 

We have not worked out theoretically how fast log z converges to \imL^^L^^oo[f {Lq+ 
l/2,t) — /(Lq- l/2,t)]. A natural guess is that the convergence is exponentially fast 
in Li,L2 in the high temperature phase of the film, while in the low temperature 
phase it follows a power law. For the thickness Lq = 8.5 we have simulated at four 
values of /? for a large range of Li = L2 up to Li = L2 = 1000. The results show 
that the convergence with Li, L2 ^ 00 is no problem in practice. Our final results 
for —A/ex at these four values of (3 are consistent but more accurate than those 
obtained in [2T] . 

Finally we have computed 9min for several thicknesses, where we have taken the 
values of Xmin from [21j. This allowed us to improve the estimate 9min = 1.31(3) 
[2T] to 9min = 1.31(2). This part of the study nicely shows that the virtues of the 
two method are complementary. 

We have not worked out theoretically how the numerical effort increases for a 
given precision with increasing thickness of the film. We also have not optimised 
the parameters of the algorithm. However it is quite clear from the simulations 
presented here that the method, using our ad hoc choice of the parameters, is 
competitive with previous proposals. 



18 



Here we have tested the method at the example of the XY universahty class. 
The application to other universality classes, like the Ising or Heisenberg universality 
class is straight forward. On the other hand, the method seems to be restricted to 
films with Dirichlet boundary conditions. 
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